clear all;
clc;

% ----------------------- Cord Geometry -----------------------------------

% spinal cord
ae = 7.44/2;         % transverse left-right axis of white matter ellipse
be = 5.5/2;          % transverse down-up axis of white matter ellipse
c_white = [0;13.61]; % center of white matter ellipse
% white matter
th_unit = 0:pi/20:2*pi;
rcord   = ae*be./sqrt((be*cos(th_unit)).^2+(ae*sin(th_unit)).^2); 
xcord   = c_white(1) + rcord.*cos(th_unit);
ycord   = c_white(2) + rcord.*sin(th_unit);
% grey matter
x1 = [0.00 ,0.62 ,1.24 ,1.87 ,2.09 ,2.24 ,1.45 ,2.03 ,2.22 ,2.03 ,1.21 ,0.62 ,0.00];
y1 = [13.17,13.49,14.57,15.65,15.67,15.49,14.04,12.49,11.99,11.49,11.43,12.00,12.73];
xgrey = [x1,-fliplr(x1)];
ygrey = [y1,fliplr(y1)];
         
% ----------------------- Importing Data ----------------------------------

relec = [0,17.38];    % 0 degrees
% relec = [0.66,17.32]; % -10 degrees
% relec = [1.29,17.15]; % -20 degrees


pol = -1;
noaxons = 200;
dcthresh = load('neg_sub_bp_dc_act_p1_d0p2_t0.txt');
[dcthresh(:,1),dc_si] = sort( dcthresh(:,1)*pol );
dc_xy = load('p1_dc_xyloc.txt');
dc_xy = dc_xy + ( c_white*ones(1,noaxons) )';
dc_xy = dc_xy(dc_si,:);

drthresh = load('neg_sub_bp_dr_act_p1_d0p2_t0.txt');
drthresh(:,1) = sort( drthresh(:,1)*pol );
dr_xy = load('p1_dr_xyloc.txt');
dr_xy = dr_xy + ( c_white*ones(1,noaxons) )';

% ----------------------- Spatial Activation ------------------------------

minth = min(dcthresh(:,1));
maxth = max(dcthresh(:,1));

frac_lv  = [0.01,0.1,0.25,0.5,0.75,0.9];
vstim_lv = minth + frac_lv*(maxth-minth);
nolev    = length(frac_lv);

figure;
for k = 1:nolev
    
    subplot(2,3,k);
    hold on;
    fill(xcord,ycord,[1,1,1],'LineWidth',3);
    fill(xgrey,ygrey,[0.5,0.5,0.5],'LineWidth',3);
    plot(relec(1),relec(2),'b.','MarkerSize',30);
    imax_dc = round( frac_lv(k)*noaxons );
    ck_dr = ( drthresh(:,1) < dcthresh(imax_dc,1) );
    plot(dc_xy(1:imax_dc,1),dc_xy(1:imax_dc,2),'g.','MarkerSize',25);
    plot(dr_xy(ck_dr,1),dr_xy(ck_dr,2),'r.','MarkerSize',25);
    ttl = [num2str(frac_lv(k)*100),' % DC'];
    title(ttl,'FontSize',30);
    xlabel('x (mm)','FontSize',20);
    ylabel('y (mm)','FontSIze',20);
    set(gca,'FontSize',20);
    hold off;
    
end
legend('White Matter','Grey Matter','Electrode',...
           'DCs Activated','DRs Activated','Location','NE');